An  AEGIS-FISST  Integrated  Detection  and 
Tracking  Approach  to  Space  Situational  Awareness 


Islam  I.  Hussein 


Kyle  J.  DeMars  Richard  S.  Erwin 

Carolin  Friih  Moriba  K.  Jah 


Worcester  Polytechnic  Institute, 
Department  of  Mechanical  Engineering, 
Worcester,  MA  01609. 

Email:  ihussein@wpi.edu 


National  Academy  of  Science, 
Air  Force  Research  Laboratory, 
Space  Vehicles  Directorate, 
3550  Aberdeen  Ave  SE, 
Kirtland  AFB,  NM  87117. 
Email:  demars.kyle@ieee.org, 
carolin@ece.unm.edu 


Air  Force  Research  Laboratory, 
Space  Vehicles  Directorate 
3550  Aberdeen  Ave  SE, 
Kirtland  AFB,  NM  87117. 


Abstract  Space  Situational  Awareness  (SSA)  is  composed  of 
three  interdependent  tasks:  discovery  of  new  objects,  tracking  of 
detected  objects,  and  characterization  of  tracked  objects.  Cur¬ 
rently  these  problems  are  treated  separately  and  independently 
of  each  other,  which  may  result  in  the  non-optimal  processing  of 
data,  with  a  corresponding  loss  of  potential  information.  Given 
the  scarcity  of  sensors  available  to  perform  SSA  missions,  it  is 
crucial  that  these  resources  be  used  as  efficiently  as  possible. 
Detection  and  classification  both  involve  estimation  over  the 
space  of  discrete  variables  (e.g.,  existence/nonexistence,  satellite 
mission  type),  whereas  tracking  involves  estimation  over  a  space 
of  continuously-valued  variables  (e.g.,  position  and  velocity).  The 
current  paper  uses  Finite  Set  Statistics  (FISST)  to  formulate  a 
hybrid  SSA  estimation  problem,  which  consists  of  simultaneously 
estimating  the  number  of  objects  and  their  tracks,  in  the 
presence  of  false  alarms,  misdetections  and  sensor  noise.  The 
main  contribution  of  the  paper  is  that,  in  order  to  reduce  the 
computational  burden  entailed  in  FISST,  we  employ  a  Gaussian 
mixture  approximation,  not  to  the  first-moment  (as  in  GM-PHD) 
of  the  full  FISST  update  equations,  but  apply  the  approximation 
directly  to  the  full  FISST  equations.  The  specific  GM  technique 
we  employ  is  the  Adaptive  Entropy-based  Gaussian-mixture 
Information  Synthesis  (AEGIS).  The  approach  is  demonstrated 
on  a  simplified  SSA  application  example. 

I.  Introduction 

About  five  hundred  thousand  objects  one  centimeter  in  size 
or  larger  populate  the  space  around  earth,  which  increases 
the  need  for  new  methods  and  techniques  that  can  contribute 
to  the  core  elements  of  Space  Situational  Awareness  (SSA): 
discovery  of  new  objects,  tracking  of  detected  objects,  and 
characterization  of  tracked  objects.  Whereas  most  SSA  lit¬ 
erature  treats  detection  [1],  characterization  [2]  and  tracking 
[3]  independently  from  each  other,  the  current  paper  seeks  a 
integrated  approach  to  all  these  problems.  In  treating  detection, 
characterization  and  tracking  separately,  valuable  information 
can  be  lost.  By  taking  into  account  the  fact  that  detection, 
characterization,  and  tracking  are  interdependent,  information 
loss  is  limited.  For  example,  object  type  may  define  some  basic 
shape  characteristics  that  affect  motion  (e.g.,  drag  coefficient 
for  low  Earth  orbit  spacecraft),  and  thus  provides  information 


about  the  object’s  dynamics  and  its  track,  and  vice  versa. 
Methods  that  explicitly  account  for  these  interdependencies 
are  therefore  at  an  advantage  over  methods  that  assume 
independence  between  the  problems. 

Any  such  integrated  approach  to  the  SSA  problem  must 
deal  with  the  hybrid  (discrete  and  continuous  variables)  in¬ 
ference  problem  that  results.  The  problems  of  detection  and 
characterization  involve  estimation  over  the  space  of  dis¬ 
crete  variables  (e.g.,  existence/non-existence,  satellite  mission 
type),  whereas  tracking  involves  estimation  over  a  space  of 
continuously-valued  variables  (e.g.,  position  and  velocity). 
Given  the  potential  coupling  between  the  dynamics  of  these 
variables  as  discussed  above,  the  SSA  problem  is  inherently  a 
hybrid  estimation  problem,  where  the  discrete  and  continuous 
variables  are  dynamically  coupled. 

Finally,  it  is  important  to  note  that  sensor  observation  time 
is  expensive,  and  only  a  few  stations  worldwide  can  acquire 
SSA  measurements.  In  particular,  sensor  resources  at  our 
disposal  are  scarce  in  comparison  to  the  size  of  the  SSA 
tasks  at  hand  (the  large  number  of  objects  to  be  catalogued 
and  the  large  search  space).  These  same  limited  sensors  can 
be  operated  in  different  sensing  modes  which  optimize  their 
detection  (for  search)  or  tracking  performance  characteristics. 
However,  wide-angle  measurements  (ideal  for  detection),  for 
example,  may  result  in  reduced  position  accuracy  and/or  loss 
of  less  reflective  (with  respect  to  sensor  detection  wavelength) 
objects.  Decisions  on  whether  to  operate  the  sensor  in  a  track 
versus  a  detect  mode  require  an  assessment  of  how  much 
information  we  expect  to  gain.  This  expected  information 
gain  will  have  two  (interdependent)  components:  one  that  is 
continuous  in  nature  and  the  other  that  is  discrete  in  nature. 
Therefore,  the  expected  information  gain  associated  with  each 
sensor  assignment  is  of  a  hybrid  nature.  The  final  requirement, 
therefore,  for  an  effective  approach  to  solving  the  SSA  prob¬ 
lem  must  provide  rigorous  mechanism  for  quantifying  hybrid 
information  gain  for  optimal  sensor  allocation. 

The  hybrid  estimation  method  we  propose  for  solving 
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the  SSA  problem  is  Finite  Set  Statistics  (FISST)  [4],  [5]. 
FISST  is  a  Bayesian  approach  that  allows  estimation  of 
hybrid  variables,  without  separating  tracking,  detection  and 
characterization.  FISST  naturally  suggests  a  way  to  generalize 
existing  measures  of  information  gain  for  purely  discrete  and 
purely  continuous  problems  to  rigorous  hybrid  measures  of 
information  gain  [5],  [6],  which  is  important  for  providing 
information-theoretic  criteria  for  sensor  allocation. 

The  greatest  challenge  in  implementing  FISST  in  real-time, 
which  is  critical  to  any  viable  SSA  solution,  is  computational 
burden.  The  first-order  Probability  Hypothesis  Density  (PHD) 
approach  has  been  proposed  as  a  computationally  tractable 
approach  to  applying  FISST  [5],  PHD  can  further  employ 
a  Gaussian  Mixture  (GM)  or  a  particle  filter  approximation 
to  further  reduce  the  computational  burden  (by  removing  the 
need  to  do  a  rigorous  discretization  of  the  state  space)  at  the 
expense  of  simplifying  assumptions.  In  this  paper,  we  develop 
a  GM  approximation  and  apply  it,  not  to  the  first-order  PHD 
approximation,  but  to  the  original  full  hybrid  propagation  and 
update  equations  derived  from  FISST.  This  eliminates  any 
information  loss  associated  with  using  the  first-order  PHD 
approximation.  The  approach  we  pursue  here  is  similar  in 
spirit  to  the  concept  of  the  “para-Gaussian”  that  was  very 
briefly  described  in  [7]. 

The  particular  GM  technique  we  use  for  the  approximation 
is  Adaptive  Entropy-based  Gaussian-mixture  Information  Syn¬ 
thesis  (AEGIS)  [8],  [9],  which  is  an  estimation  approach  for 
non-linear  continuous  dynamical  systems.  AEGIS  implements 
a  GM  model  representation  of  the  probability  density  function 
(pdf)  that  is  adapted  online  via  splitting  of  the  GM  components 
whenever  an  entropy-based  detection  of  nonlinearity-induced 
distortions  of  the  Gaussian  components  is  triggered  during 
the  forward  propagation  of  the  pdf.  In  doing  so,  the  GM 
approximation  adaptively  includes  additional  components  as 
nonlinearity  is  encountered  and  can  therefore  be  used  to  more 
accurately  approximate  the  pdf. 

To  summarize,  the  main  contributions  of  this  paper  are 
(1)  the  use  of  the  GM  AEGIS  approach  to  approximate  the 
complete,  un-approximated  FISST  propagation  and  update 
equations,  and  (2)  the  application  of  this  approach  to  the 
problem  of  detection  and  tracking  of  space  objects.  The 
latter  is  demonstrated  via  a  simple  “0-1”  object  example.  The 
remainder  of  the  paper  is  organized  as  follows.  We  first  briefly 
summarize  the  theory  behind  FISST  in  Section  II.  Next  we 
briefly  summarize  AEGIS  in  Section  III.  In  Section  III-A  the 
general  AEGIS-FISST  approach  is  summarized.  In  Section  IV, 
the  method  is  applied  to  a  planar  single  object  SSA  detection 
incorporating  the  possibility  of  false  alarms,  misdetections, 
and  sensor  tracking  noise.  In  Section  V  we  present  numerical 
simulation  results  for  the  simple  SSA  problem,  demonstrating 
the  capabilities  and  performance  properties  of  the  method.  We 
conclude  with  future  research  directions  in  Section  VI. 

II.  A  Brief  Introduction  to  FISST 

As  opposed  to  purely  discrete  or  purely  continuous  Bayesian 
inference,  FISST  makes  use  of  set-valued  random  variables. 


For  example,  if  we  let  W  be  the  set  of  all  possible  space 
object  types,  then  Xd  £  W  is  the  discrete  component  of  the 
state  that  describes  a  space  object’s  type  and  x  g  l6  is  the 
continuous  component  of  the  state  that  specifies  the  position 
and  velocity  of  the  object,  then  the  set-valued  random  variable 
corresponding  to  the  hybrid  system  state  is  X  =  {x,j,  x}. 
In  a  detection  and  tracking  problem,  the  set-valued  hybrid 
system  state  X  =  (n.  X),  where  n  is  the  discrete  component 
that  describes  the  number  of  objects  in  the  search  space  and 
X  =  [x^  xj  ...  xj]  £  M6"  describes  the  positions  and 
velocities  of  these  objects.  Notice  here  the  explicit  dependence 
of  the  dimension  of  the  continuous  state  space  K.6rl  on  the 
discrete  component  n  of  the  state.  For  brevity,  we  simply  write 
X  =  {xi,  X2, . . . ,  x„  }  in  the  detection  and  tracking  problem. 

In  this  paper  we  only  look  at  the  integrated  detection  and 
tracking  problem.  Hence,  the  state  X  =  {xi,X2, . . . ,  xn} 
gives  the  number  of  objects  n  and  where  x,  describes  the 
position  and  velocity  of  object  i.  We  will  assume  a  planar 
dynamic  model.  In  other  words,  we  have  x,  £  K4.  Similar 
definitions  apply  to  the  set-valued  measurement  variable  Z  = 
{zi,  Z2, . . . ,  zmj  with  to  sensor  returns  and  where  z £  Rm, 
i  =  1, . . . ,  to,  is  the  value  of  each  return. 

Bayes’  law  takes  up  exactly  the  same  form  in  the  hybrid 
FISST  approach  as  it  does  in  purely  continuous  and  purely 
discrete  problems: 

/fc+1|fe(X|Z(fc))  =  f  fk+Mk(X,X')  ■  fklk(X'\zW)6X' 

(  t  V\  f7(k+l)\  fk+l(Zk+l  |^D/fc  +  l|fc  (-Xj-Z^  /1\ 

h+i\k+i{x\z  )-  fk+1(zk+1\z(^)  >  W 

where  Z (k+1)  =  {Z\,...,  Zk+\]  is  the  time  sequence  of 
measurement  sets  up  to  time  k  +  1,  and  where  /fc+i|fc(X|X') 
is  the  multi-target  Markov  transitional  density.  The  function 
fk+i(Z\X)  is  the  multi-target  likelihood  function  that  de¬ 
scribes  the  likelihood  of  getting  a  measurement  Z^+i  given 
the  state  X^+i-  The  first  of  Eq.  (1)  is  the  prediction  step  and 
the  second  is  the  update  or  corrector  step.  The  normalization 
factor  in  the  update  step  is  given  by 

fk+1(Zk+1\zW)  =  j  fk+1(Zk+1\X)fk+1{k(X\Z^)SX.  (2) 

Notice  that  the  integrals  are  set  integrals,  denoted  by  the 
variable  of  integration  SX  vice  dX  for  a  standard  integral. 
For  multi-target  detection  and  tracking,  the  set  integral  of  a 
scalar-valued  set  function  g(X)  is  defined  to  be  the  integral 
of  g  over  the  continuous  component,  summed  over  all  possible 
discrete  values  [4],  [5]: 

J  g(X)SX  =  g(X  =  0) 

+  X!"T  /  ff({xi,...,Xi})dxi  •••dxi(3) 

i= 1 

where  the  factorial  coefficient  takes  into  account  all  the 
different  possible  orderings  of  X  as  evaluated  in  the  function 
g.  The  need  for  this  factor  accounts  for  the  data-association 
problem  [5]. 
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III.  A  Brief  Introduction  to  AEGIS 

Most  implementations  of  Gaussian  mixture  propagation 
algorithms  assume  that  the  component  weights  remaining 
constant  over  the  duration  of  the  propagation  span,  and  so 
do  not  incorporate  a  method  for  online  refinement  of  the 
Gaussian  mixture.  Methods,  such  as  the  Adaptive  Gaussian 
Mixture  provide  a  method  for  adaptive  weight  variation  with¬ 
out  a  mechanism  to  vary  (coarsen  or  refine)  the  number 
of  Gaussian  components  [10].  The  AEGIS  filter  approaches 
the  problem  of  adapting  the  weights  of  the  GM  pdf  by 
monitoring  nonlinearity-induced  distortions  of  the  Gaussian 
components  during  the  propagation  of  the  pdf  and  using  a 
splitting  algorithm  to  increase  the  accuracy  of  linearization, 
thereby  allowing  the  filter  to  modify  the  GM  components 
(weight  value  and  cardinality)  in  such  a  way  so  as  to  avoid 
significant  linearization  errors. 

1 )  Dynamical  System  Model:  Many  systems  of  interest  fall 
under  the  broad  classification  of  nonlinear  systems.  An  estima¬ 
tion  algorithm  which  exploits  at  least  some  characteristics  of 
the  nonlinearities  is  preferable  to  approximating  the  problem 
as  that  of  a  linear  one.  Consider  the  nonlinear  dynamical 
system  governed  by  the  differential  equation 

x(t)  =  f(x(f),  t) ,  x(f0)  =  X0  ,  (4) 

where  x(i)  is  the  non-set  valued,  continuous  state  of  the 
system,  f(-)  represents  the  sufficiently  differentiable  nonlinear 
dynamics  of  the  system,  and  xo  is  the  initial  condition.  The 
initial  condition  is  assumed  to  be  random  with  pdf  /0|0(x). 

2)  Detecting  Nonlinearity  Induced  Distortions  to  the  GM 
during  Propagation:  An  integral  aspect  of  the  AEGIS  propa¬ 
gation  scheme  is  the  detection  of  nonlinearity  during  the  prop¬ 
agation  of  uncertainty.  The  method  employed  in  the  AEGIS 
approach  is  based  on  a  property  derived  from  the  differential 
entropy  for  linearized  dynamical  systems  that  allows  for 
the  detection  of  nonlinearity  without  employing  linearization- 
based  methods.  It  can  be  shown  that  the  differential  entropy 
for  a  linearized,  Gaussian  system  evolves  as  [11],  [12] 

H(x)  =  trace {F(x(f),f)}  ,  (5) 

where  F(x(t),t)  is  the  dynamics  Jacobian  matrix.  The  value 
of  the  entropy  for  a  linearized  system  can  be  determined 
by  numerically  integrating  Eq.  (5)  with  an  appropriate  initial 
condition  and  requiring  only  the  evaluation  of  the  trace  of 
the  linearized  dynamics  Jacobian.  In  parallel,  a  nonlinear 
implementation  of  the  integration  of  the  covariance  matrix 
(e.g.  the  Unscented  Kalman  Filter  (UKF)  [13])  is  considered, 
which  allows  a  nonlinear  determination  of  the  differential 
entropy.  Any  deviation  between  these  values  of  entropy  there¬ 
fore  serves  as  an  indication  that  nonlinearity  is  impacting  the 
solution.  When  the  difference  between  the  linearized-predicted 
entropy  and  the  nonlinear  computation  of  the  entropy  exceeds 
a  user-defined  threshold,  nonlinearity  has  been  detected  in  the 
propagation  of  the  uncertainty,  the  propagation  is  halted,  a 
splitting  algorithm  is  applied  to  the  Gaussian  distribution,  and 
propagation  resumes  with  the  adapted  distribution. 


3)  Propagation:  Consider  the  time -propagation  of  the  pdf 
and  consider  the  time  interval  t  £  [ifc,  tfc+i]-  It  is  desired  to 
approximate  the  conditional  pdf  at  time  tk  via 

N 

A-|fc(x|^fe)  =  '^2wiNi(x;p.i,Pi),  (6) 

2—1 

where  TV^x;  //, ,  P.j)  is  a  Gaussian  distribution  in  x  with  mean 
jli  and  covariance  P,,  and  where  w.L  is  the  weight  of  the  ith 
GM  component  with  =  1-  The  propagated  pdf  is 

then  given  by 

N 

fk+i\k{*\Zk)  =  (7) 

2=1 

where  uii,  fi,  and  P,,  i  =  1  are  the  weight,  mean 

and  covariance  of  the  of  the  ith  propagated  GM  component. 
It  should  be  noted  that  due  to  the  splitting  algorithm  described 
above,  the  number  of  components  in  fk+i\k(x\Zk),  given  by 
N,  may  now  be  different  than  the  number  of  components  in 
A|fc(x|Zfc),  given  by  N. 

To  propagate  the  pdf  forward  in  time,  a  nonlinear  algorithm 
such  as  the  UKF  is  applied  to  each  component  of  the  GM  pdf. 
Simultaneously,  the  linearized  differential  entropy  is  propa¬ 
gated  as  described  previously.  Each  component  is  monitored 
for  deviations  in  the  nonlinear  and  linear  predictions  of  the 
entropy.  Once  nonlinearity  is  detected  on  a  component,  the 
propagation  is  halted  at  time  is,  where  tk  <  ts  <  tk+i-  If 
ts  /  tk+i,  then  a  splitting  step  is  performed  on  the  component 
for  which  nonlinearity  was  detected.1  That  is,  if  nonlinearity 
was  detected  in  the  jth  component  at  time  ts,  then  the  jth 
component  is  replaced  at  time  ts  by 

G 

utjNj(x;  «  y ^wrNr(x-  jttr,Pr)  (8) 

r—  1 

where  ts  has  been  omitted  for  brevity,  and  the  replacement 
component  weights,  means,  and  covariances  are  computed 
using  the  Gaussian  splitting  algorithm  discussed  in  [11].  After 
the  splitting  step  has  been  performed,  return  to  Eq.  (7)  with 
tk  =  ts  and  N  ■£-  N  +  G  —  1  components  in  the  Gaussian 
mixture,  then  continue  until  ts  =  tk+i  is  reached.  Once 
ts  =  tk+ i,  the  propagation  step  has  been  completed  with  N 
components  having  weights  uii,  means  fi,,  and  covariances 
Pi,  which  allows  the  a  priori  GM  pdf  to  be  evaluated  via 
Eq.  (6). 

A.  AEGIS  based  Approximation  of  FISST:  AEGIS  FISST 

The  AEGIS-FISST  approach  relies  on  the  observation  that 
if  the  (conventional  non-hybrid)  prior  fk\k{~x\Z^)  is  a  GM, 
then  the  posterior  multi-target  pdf  will  also  be  a  sum  of 
Gaussian  mixtures  from  which  one  can  extract  the  (conven¬ 
tional)  posterior  pdf  ,fk+i\k+i(x\Zf'k+1'>)  (also  a  GM)  under 
the  hypothesis  that  an  object  exists.  This  observation  is  what 

1  Without  loss  of  generality,  nonlinearity  is  assumed  to  be  detected  on  only 
one  component.  If  more  components  detect  nonlinearity,  the  same  process  is 
applied  to  each  component  individually. 
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distinguishes  the  proposed  method  from  a  GM  implementation 
of  the  PHD  approach  [5].  PHD  relies  on  approximating  the 
multi-target  posterior  pdf  by  its  first-order  moment  PHD. 
After  that,  a  GM  is  used  to  approximate  the  PHD.  In  our 
approach,  we  circumvent  the  first  step  of  approximating  the 
multi-target  pdf  by  its  PHD  and  realize  that  the  Gaussian 
mixture  property  is  preserved  in  the  update  step.  Thus,  one 
need  not  approximate  the  multi-target  pdf  by  its  first-order 
PHD  and  the  only  approximation  one  performs  is  that  of  using 
a  GM  to  model  the  prior  pdf. 

IV.  A  Single-Object  Illustrative  Example 
A.  The  Basic  FISST  Propagation  and  Update  Equations 

For  the  sake  of  brevity  and  ease  of  exposition,  we  consider 
a  problem  in  which  there  can  exist  at  most  one  object  in  the 
search  space.  The  sensors  are  assumed  imperfect  with  a  non¬ 
zero  probability  of  false  alarm  due  to  clutter  (small  objects  that 
are  of  no  interest,  thruster  exhaust  plumes,  environmentally- 
induced  sensor  responses,  etc),  and  with  a  non-unity  proba¬ 
bility  of  detection  when  the  object  is  within  a  pre-specified 
field-of-view,  and  a  zero  probability  of  detection  when  the 
object  is  outside  this  field-of-view.  In  tracking  an  object,  the 
sensor  adds  some  tracking  noise.  Without  loss  of  generality, 
the  sensor  is  assumed  to  measure  the  position  and  velocity  of 
an  object  directly. 

For  the  single  object  case  with  at  most  a  single  source  of 
clutter,  one  can  use  a  Bernoulli  distribution  [14]  for  the  various 
multi-target  densities  as  we  will  show  below.  Extension  to  the 
arbitrary  number  of  objects  can  be  addressed  using  a  Poisson 
distribution  model  for  the  number  of  objects  as  well  as  clutter 
sources. 

Under  these  assumptions  and  given  measurements  Zf  k)  up 
to  time  step  k,  the  multi  target  prior  density  function  has  the 
form 

[  1  —Pk  if  X  =  0 

fMk(X\zW)  =  {  Pk  ■  fk\k(x\Z^)  ifX  =  {x}(9) 

[  0  if  |.Xj  >  2 

where  pk  is  the  prior  probability  that  the  object  exists  at 
time  k  (we  will  henceforth  drop  the  subscript  k)  with  prior 
conventional  pdf  given  by  fk\k(x\Z^). 

Remark.  In  composing  a  multi-target  pdf,  say  for  the  the 
case  that  X  =  {x}  (i.e.,  the  hypothesis  that  an  object 
exists  in  the  search  space  with  a  state  x),  one  considers 
fk\k(X  =  {x}|  ZW)  as  the  “probability”  of  an  object  ex¬ 
isting  with  true  probability  p  and  having  the  object’s  con¬ 
tinuous  state  having  (conventional)  pdf  fk\k(x\Z^).  Hence, 
fk\k(X\ZW)  =  p  ■  fk\k(x\Z^)  as  shown  in  the  second 
line  in  Eq.  (9).  The  “and”  operation  was  translated  into  a 
product  of  the  discrete  probability  p  and  the  continuous  pdf 
/fc|fc(x|Z(fe)).  Likewise,  the  probability  that  X  =  0  is  thus 
simply  1  —p.  One  can  check  that  J  fk\k{X\Z^)8X  =  1  and, 
hence,  fk\k{X\Z^k>)  is  a  valid  multi-target  pdf. 

Since  we  make  the  simplifying  assumption  that  at  most 
a  single  object  exists  in  the  search  space,  the  multi  target 


transitional  Markov  density  function  is  given  by 

fk+i\k(X\</))  =  |  Q  if  \X\  >  1 

(  1  -ps  if  AG  =  0 

/fe+i|fc(X|{x'})  =  <  Ps  ■  /fc+i|fc(x|x')  if  x  =  {x}  (10) 
[  0  if  |.Xj  >  2 

where  ps  is  the  probability  that  the  object  servives  in  the 
search  space  from  time  step  k  to  time  step  k  +  1  and  where 
A-+i|fc(x|x')  is  the  conventional  Markov  transitional  density 
function  that  describes  the  likelihood  of  transitioning  to  the 
state  x  from  the  state  x'  under  the  hypothesis  that  an  object 
with  state  x'  exists.  One  can  add  a  birth  model,  but  we  omit 
this  here  for  the  sake  of  transparency  of  exposition. 

Applying  the  basic  definition  of  a  set  integral,  for  the 
predicted  multi  target  density  function,  we  have 

fk+i\k(X  =  0|Z(fc)) 

=  J  fk+i\k(X  =  <b\X')fklk(X'\zW)8X' 

=  fk+i\k(m  ■  fk\k(9\Zw) 

+  J  /fc+i|fc(0|{x'}  •  /*|fc({x,}|Z(fc))dx' 

=  (1  ~P)+P  J (1  ~  Ps)fk\k{x'\Z{k))dx' 

=  (f  ~P)  +P(f  ~Ps)  (11) 


and  (we  have  omitted  the  calculation  for  the  second  expression 
for  the  sake  of  brevity) 


fk+i \k(X  =  {x}|  ZW)  =  pspfk+1\k(x\Z^).  (12) 


One  can  check  that  f  fk+i\k(X)SX  =  1. 

Let  the  probability  of  detection  be  pr>  and  the  probability 
of  false  alarm  be  pp.  These  two  probabilities  are  generally 
a  function  of  the  location  of  the  object  (whether  it  is  in  the 
field-of-view  or  not).  Then  the  multi  target  likelihood  function 
is  given  by 


fk+l{Zk+l  =  0|0) 

fk+ l(Zk+l  =  {z}|0) 
fk+i(Zk+1  =  0|{x}) 
fk+i{Zk+1  =  {z}|{x}) 

fk+i{Zk+ 1  =  {zi,z2}|{x}) 


1  -  Pf 

PFg(z) 

(1  -  pF){l  -  Pd) 

Pf{  1  ~PD)g{?) 

+Pd{  1  -PF)/fc+i(z|x) 
PfPd  [f/(zl )  A+l  (z2  |x) 
+ff(z2)/fc+l(zi|x)]  (13) 


where  g( z)  is  the  spatial  likelihood  clutter  distribution  function 
that  a  clutter  generates  the  measurement  z. 

Remark.  In  order  to  understand  where  the  components  of 
the  multi-target  likelihood  function  originate  from,  consider 
for  example  the  last  component  fk+i (Zk+i  =  {zi,  Z2 } | {x}). 
Given  that  an  object  exists  in  the  search  space,  this  is  the 
“probability”  that  two  measurements  are  registered.  The  only 
two  possibilities  are  that  zi  was  generated  by  a  clutter  source 
and  z2  by  an  object,  or  vice-versa.  This  gives  the  sum 
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in  the  bracketed  term.  For  that  scenario  (i.e.,  two  sensor 
returns)  to  occur  a  detection  of  the  object  (hence,  the  pn 
coefficient  outside  the  bracketed  term)  as  well  as  a  false 
alarm  triggered  by  clutter  (hence,  the  p p  coefficient  outside 
the  bracketed  term)  must  take  place.  One  can  check  that 
J  fk+i(Zk+i\$)5Zk+i  =  f  fk+i(Zk+1\{x})SZk+1  =  1. 

The  multi  target  Bayes  factor  in  the  denominator  of  the 
corrector  step  (the  second  of  Eq.  (1)  is  given  by 

fk+i{Zk+i  =  0|  Z^) 

=  I  fk+1(Z  =  ®\X)fk+llk(X\zW)6X 

=  (!  ~Pf)  ((1  -p)  +p(l  ~Ps))  +  (1  ~Pf){  1  ~  Pd)PsP, 

(14) 


and 

/fc+i|fc+i(^  =  {x}\Zk+i  =  {z})  = 

(pD(l  -pF)/fc+i(z|x)  +pF(  1  -  Pd)s(z)) PsPfk+Hki^Z^) 

fk+i(Zk+i  =  M) 

(20) 

If  Zk+ 1  =  {zi,z2}  we  have: 

fk+i \k+i(X  =  $\Zk+i  =  {Zl,z2})  =  0.  (21) 

Remark.  Equation  (21)  makes  perfect  sense  since  getting  two 
returns  (under  the  assumption  that  at  most  a  single  object 
exists)  ensures  that  at  least  one  of  them  was  generated  by  an 
existing  object  (and,  hence,  that  the  probability  that  no  object 
exists  is  zero),  and 


fk+1(Zk+1  =  {z}\Z^)  =  J  fk+1({z}\X)fk+llk(X\zW)6X 

=  [: PF  ((1  -p)+p{  1  -  Pa))  +PsPPf{  1  -pz>)]s(z) 

+  PsPPd(1  -pF)fk+i(z),  (15) 

and 

fk+l(Zk+l  =  {zuz2}\Z^) 

=  J  fk+1(Z  =  {z1,z2}\X)fk+1{k(X\Z^)8X 

=  PsPPfPd  (g( Zi)/(z2)  +  s(z2)/(zi)) ,  (16) 

where  fk+  i(z)  is  the  conventional  Bayes  factor  given 
that  an  object  exists.  Again,  one  can  check  that 

f  fk+1(Zk+1\ZW)8Zk+1  =  l. 

Putting  all  the  above  together  for  the  multi  target  corrector 
step,  we  get  the  following. 

If  Zk+ 1  =  0  (i.e.,  no  sensor  return  at  time  step  k  +  1) 
the  posterior  multi-target  pdf  component  corresponding  to  the 
hypothesis  that  no  object  exists  is  given  by: 

/fc+i|fc+i(*  =  0|^fc+i  =  0) 

=  fk+i(Zk+1  =  0|Xfc+1  =  <D)fk+i\k(X  =  0|  Z^) 

fk+l(Zk+i  =  0) 

_  (1  -  Pf)  ((1  -  p)  +  p(l  -  Ps))  n 

fk+i(Zk+i  =  0) 


fk+i\k+i(X  —  {x}|^fc+i  —  {zi,z2})  — 

pFpg(g(zi)/fc+i(z2|x)  +g(z2)/fc+i(zi|x))psp/fc+1|fc(x|Zlfcl) 
/fe+l(^fe+l  =  {zi,z2}) 

(22) 

This  second  component  will  in  fact  integrate  to  one.  The 
integral  of  fk+1\k+1{X  =  {x.}\Zk+1  =  {zi,z2})  gives  the 
posterior  discrete  probability  that  an  object  exists.  Hence,  that 
quantity  will  have  to  integrate  to  1. 

B.  The  AEGIS  FISST  Approximation  for  the  Single  Object 
Problem 

We  will  assume  standard  (unperturbed)  Keplerian  object 
dynamics  which  can  be  modeled  in  the  form  (4).  Substituting 
(6)  and  (7)  into  Eq.  (17)-(22)  to  give  an  updated  set  of  posterior 
Gaussian  sum  parameters  at  time  fc+1.  The  posterior  Gaussian 
sum  has  the  form 

N 

fk+1\k+1(x\Zk+1)  ='^2wiNi(x-,j. hi,  Pi) 

i= 1 

_  fk+\\k+l{X  ~  Ml-Zfc+l) 

P(tk+ 1)  ’  (  } 

with  Eili  Wi  =  i  and  where  p(tk+i)  is  the  discrete  posterior 
probability  that  an  object  exists  and  is  given  by  [5] 

p(tk+i)=  J  fk+i\k+i{X  =  {x}|Zfc+i)dx.  (24) 


Likewise  for  the  case  that  an  object  exists  with  state  x,  we 
obtain 


fk+l\k+l(X  —  {x}|Zfc+i  —  0) 

_  (1  ~P£>)(1  -  PF)PsPfk+l\k(^\Z('k)) 

fk+l(Zk+l  =  0) 


(18) 


Omitting  the  detailed  derivation,  if  Zk+\  =  {z}  we  have: 


fk+i\k+i(X  —  0|^fc+i  —  {z}) 

=  pf  ((l  -p)+p( l  -Ps))g{ z) 

fk+l(Zk+l  =  {z})) 


(19) 


Equation  (23)  can  be  used  to  extract  the  components  of  the 
posterior  GM  under  the  hypothesis  that  an  object  exists. 

If  Zk+ 1  =  0,  since  there  is  no  return,  the  “updated” 
components  are  simply  the  same  as  those  of  the  propagated 
GMs.  This  can  readily  be  verified  by  substituting  Zk+ 1  =  0 
into  the  update  equations  from  the  previous  section.  Hence, 
the  updated  GMs  are  given  by: 

N  =  N, 

Wi  =  Wi, 

pi  —  Pi 

P*  =  Pi,  (25) 
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where 


i  =  1, . . .  ,N. 

If  Zk+ 1  =  {z}  (i.e.,  the  sensor  gets  one  return),  the  resulting 
posterior  GM  can  be  discretized  down  into  two  separate  sums: 
(1)  an  update  of  the  posterior  under  the  hypothesis  that  the 
measurement  came  from  an  existing  object,  and  (2)  an  update 
under  the  hypothesis  that  the  measurement  came  from  a  clutter 
and,  in  which  case,  the  propagated  GM  is  left  unchanged  for 
this  case.  Total  number  of  components  would  then  be  IV  = 
2 N.  Substituting  Zk+i  =  {z}  and  Eq.  (23)  into  the  update 
equations  for  the  single-object  problem,  one  can  check  that 
the  posterior  components  of  the  GM  are  given  by  (we  split 
the  components  into  two  groups  for  the  two  cases  pointed  out 
earlier  in  the  paragraph): 

N  =  2N 

Pd{  1  -  pF)pspNi  fz;  H/ij,  R  +  HP;HT)  u>i 
2  di(z) 

-2  _  PF(l-PD)PsPg(z)Wi 

i  di(z) 

p\  =  (P"1  +  HTR_1H)_1(P“1/ij  +  HTR-Y) 

Pi  =  Pi 

Pi  =  (p-'  +  HTR-'h)”1 
P?  =  Pi, 

i  =  1 , . . . ,  N,  where 

dl(z)  =  Pd{\  -  PF)pspfk+l(z)  +PfO-  ~  PD)pspg(z), 


d2(z1,z2)  =  PfPd  (g(zi)pspfk+i(z2)  +  g(z2)pspfk+1(z1)) , 

and  where  fk+i  (z7).  j  =  1,2  for  the  Gaussian  sum  approxi¬ 
mation  can  be  shown  to  be 

N 

fk+1  {Zj )  =  Y,  (z j  ;  H Ai ,  R  +  HP  Ht  ) .  (28) 

i=l 

The  above  posterior  Gaussian  sum  parameters  become  the 
prior  at  time  step  k+ 1.  The  above  GM  parameters  can  be  used 
to  compute  the  posterior  multi-target  pdf  fk+i\k+i(X\Zk+i). 
One  can  then  employ  the  Marginal  Multitarget  (MaM)  esti¬ 
mator  or  the  Joint  Multitarget  estimator  to  solve  for  the  state 
estimate  xj.+1  if  an  object  is  determined  to  be  existing  [5].  For 
this  simple  single-object  case,  both  methods  are  equivalent 
and  the  state  estimate  Xfc+i  is  simply  the  the  maximum  of 
fk+i\k+i(X\Zk+i)- 

In  general  detection  and  estimation  problems,  the  first  step 
in  the  MaM  estimator  is  to  compute  the  maximum  a  posteriori 
(MAP)  estimate  of  the  number  of  objects,  which  is  given  by: 

n  =  round(argsup/fc|fc(n|Z(fe))),  (29) 

n 

where  /*.!/.  (n|z?(fc))  is  the  cardinality  distribution  [5].  Given 
(26)  n,  the  second  step  is  to  solve  for  the  maximum  a  posteriori 
estimate  of  the  state 

xMaM  =  arg  sup  /fc|fc({xi,...,x*|Z(fc)}).  (30) 


and  where 


N 

fk+i( z)  =  ^  T25jJVi(z;  H£i,  R  +  HP,HT). 

i= 1 


In  the  above  equations  H  is  the  linearization  of  h{x)  evaluated 
at  the  propagated  value  of  the  state  estimate  x,  . 

If  Zk+ 1  =  {zi ,  z2 }  (i.e.,  the  sensor  receives  two  returns), 
the  resulting  posterior  Gaussian  sum  can  be  broken  down  into 
two  separate  sums:  (1)  an  update  of  the  posterior  under  the 
hypothesis  that  the  first  measurement  came  from  a  clutter  and 
the  second  from  an  existing  object,  and  (2)  an  update  of  the 
posterior  under  the  hypothesis  that  the  second  measurement 
came  from  a  clutter  and  the  first  from  an  existing  object.  The 
total  number  of  components  in  the  posterior  GM  would  then 
be  IV  =  2 N  with  components 


N 


w] 


Y2 


Mi 

Pi 

P  l 


2  N 

PFPDPsPg{zi)Ni  (z2;  H£j,  R  +  HP!,HT^  wt 
d2(  zi,z2) 

PFPDPsPg{z2)Ni  (zx-Ujii,  R  +  HP4Ht)  Wi 

d2(  zi,z2) 

(P-1  +  HTR_1H)_1(P~1/ij  +  HtR”1z2) 
(P"1  +  tf’R-^)-1^1/^  +  hTR-^O 
p2  =  ^p-1  +  HTR-1^  1 ,  (27) 


The  MaM  estimator  is  Bayes-optimal  but  is  not  known  if  it  is 
statistically  consistent  (i.e.,  whether  it  converges  to  the  correct 
state  in  the  static  case)  [4], 

For  the  (at  most)  single  object  object  problem,  n  =  1  if 
the  posterior  discrete  probability  p(tk+  i)  is  larger  than  0.5.  If 
p(tk+ 1)  >  0.5,  the  state  estimate  corresponds  to  the  maximizer 
of  the  posterior  multi-target  pdf  fk+i\k+i(X  =  {x}|Zfc+i). 

V.  Numerical  Simulation  Results 

In  this  section  we  demonstrate  the  single-object  approach, 
outlined  in  the  previous  section,  on  a  simple  SSA  problem. 
We  assume  that  there  can  be  at  most  a  single  object.  If  the 
orbiting  object  exists,  it  satisfies  unperturbed  two-body  planar 
Keplerian  dynamics  [15].  We  will  assume  that  if  an  object 
indeed  exists,  its  probability  of  survival  is  ps  =  1. 

We  set  the  problem  up  such  that  an  object  initially  exists 
with  probability  one  and  set  the  prior  assumed  by  AEGIS- 
FISST  to  be  p0  =  1.  The  object  true  initial  state  is  such 
that  its  apoapsis  is  42,100  km  and  the  eccentricity  is  0.2 
(specifying  these  determines  the  initial  position  and  velocity 
of  the  object  [15]).  The  prior  distribution  fed  into  the  AEGIS- 
FISST  algorithm  was  assumed  to  have  a  mean  identical  to  the 
true  object  location,  but  with  a  position  variance  of  1  km2 
and  a  velocity  variance  of  1  m2/s2.  Hence,  the  initial  prior  is 
a  single-component  Gaussian  mixture  distribution  as  required 
by  AEGIS-FISST. 
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The  clutter  distribution  g( z)  at  time  that  generates  false 
alarms  is  assumed  to  be  Gaussian  centered  around  the  current 
location  and  velocity  of  a  small,  insignificant  object  such  that: 
g( z)  ~  JV(z;xc,  Rc),  where  xc  is  the  state  of  the  clutter 
source  and  is  assumed  to  satisfy  (unperturbed)  Keplerian 
dynamics  itself  with  apoapsis  of 42,000  km  and  an  eccentricity 
of  0.6.  The  clutter  measurement  error  covariance  is  given  by: 


Re 


0.5  0  0  0 

0  0.5  0  0 

0  0  lO"3  0 

000  10"3 


Modeling  clutter-generation  as  the  result  of  a  single  clutter 
object  source  in  orbit  around  Earth  is  an  assumption  that  we 
introduce  for  the  sake  of  simplicity.  A  more  realistic  Poisson 
model  generating  multiple  clutter  points  can  be  introduced  for 
increased  modeling  fidelity  at  the  cost  of  increased  simulation 
numerical  burden. 

The  sensor  is  assumed  to  be  located  on  the  surface  of  the 
Earth,  and  is  hence  rotating  with  respect  to  the  inertial  frame. 
The  sensor  field-of-view  is  assumed  to  be  ±20  degrees  from 
the  local  vertical.  If  the  object  is  within  the  field-of-view,  its 
probability  of  detection  po  =  0.6,  and  zero  otherwise.  If  the 
clutter  source  is  within  the  field-of-view,  the  probability  of 
false  alarm  is  pp  =  0.3,  and  zero  otherwise.  The  sensor 
position  tracking  error  is  assumed  to  be  zero-mean  with 
covariance  0.500  km2  and  we  assume  that  the  sensor  measures 
the  cartesian  components  of  velocity  corrupted  by  zero-mean 
noise  with  a  variance  of  1  m2/s2.  Hence,  die  output  matrix  H 
is  simply  the  identity  matrix  and 
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At  the  end  of  each  time  step,  using  the  posterior  GM 
components,  we  compute  the  posterior  multi-target  pdf 
fk+i\k+i(X \Zk+i)  and  use  Matlab’s  unconstrained  maximiza¬ 
tion  function  fminunc  to  find  the  MaM  state  estimate. 

The  simulation  was  run  for  twelve  orbital  periods  (approxi¬ 
mately  9  days)  with  a  sensor  sampling  rate  of  a  single  measure¬ 
ment  every  200  seconds.  Figure  1  shows  the  relative  position 
error  between  the  estimated  and  true  object  position.  In  the 
figure,  a  red  star  indicates  the  times  at  which  a  measurement 
of  the  object  alone  has  been  taken.  A  blue  star  indicates 
that  a  clutter-generated  measurement  has  been  recorded  (of 
course,  the  algorithm  does  not  know  whether  the  measurement 
was  generated  by  clutter  or  by  the  real  object).  A  black  star 
indicates  that  both  the  clutter  as  well  as  the  object  have  both 
been  recorded  by  the  sensor.  Times  for  which  no  stars  appear 
are  times  during  which  no  measurements  were  returned  either 
the  clutter  source  as  well  as  the  object  were  out  of  the  field-of- 
view  or  the  object  and  clutter  source  were  in  view  but  have  not 
have  generated  any  returns  (i.e.,  misdetections  without  false 
alarms).  The  velocity  error  is  shown  in  Figure  2 

Figure  3  shows  the  posterior  probability  p(tjt+j)  that  an 


object  exists  as  a  function  of  time.  Figure  4  shows  the  number 
of  GM  components  generated  by  AEGIS.  The  number  of 
components  tended  to  be  large  in  the  initial  stages  of  the 
simulation,  but  as  more  measurements  of  the  object  are  taken, 
a  fewer  number  of  GM  components  is  generated.  Even  when 
the  components  are  larger  in  number,  the  simulation  takes, 
typically,  less  than  1  second  per  simulation  time  step  (the 
entire  code  was  written  in  Matlab  and  run  on  a  2.4GHz 
computer). 

Figure  5  shows  the  hybrid  information  divergence  gained 
at  each  time  step.  We  used  a  FISST  version  of  the  gamma 
information  divergence  [16]  to  evaluate  the  gain.  How  to 
compute  the  gamma  divergence  in  FISST  is  discussed  in 
[6].  As  we  can  see  from  the  figure,  a  very  large  amount 
of  information  is  gained  the  first  time  we  detect  the  object. 
As  long  as  the  object  is  in  the  field-of-view,  information  is 
gamed  as  tracking  measurements  are  made.  Hence,  the  jump  in 
information  divergence  is  due  to  the  (re)detection  of  the  object 
(i.e.,  a  gain  in  information  related  to  the  discrete  component 
of  the  problem)  and  the  subsequent  gain  in  information  is  due 
to  the  tracking  of  the  object’s  state  (i.e.,  a  gain  in  information 
related  to  the  continuous  component  of  the  problem).  Once 
the  object  is  no  longer  in  the  field-of-view,  the  information 
gain  is  zero. 


Fig.  1 .  Magnitude  of  the  true  relative  position  error. 

VI.  Conclusion 

In  this  paper  we  used  FISST  to  solve  a  simple  SSA  search, 
detect  and  track  problem.  A  hybrid  estimation  technique  like 
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Fig.  2.  Magnitude  of  the  true  relative  velocity  error. 
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Fig.  3.  Probability  of  object  existence. 
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Fig.  5.  Information  divergence. 


FISST  is  important  to  a  problem  like  SSA  since  information 
loss  due  to  the  artificial  separation  of  the  detection  and  tracking 
problems  is  eliminated.  However,  FISST  is  computationally 
expensive  unless  some  approximation  technique  is  employed. 
Unlike  the  GM-PHD  approach,  the  multi-taiget  pdf  is  first 
approximated  by  its  first  moment  PHD  which,  in  turn,  is 
further  approximated  using  a  GM  model.  In  the  approach 
presented  in  this  paper,  we  do  not  rely  on  using  the  PHD  as  a 
first  moment  approximation  of  the  multi-target  pdf.  Instead, 
it  is  shown  that  a  prior  pdf  that  is  represented  by  a  GM 
implies  that  the  posterior  pdf  can  also  be  represented  by  a 
GM.  This  should  reduce  the  computational  burden  by  using  a 
GM  approximation  to  FISST,  but  without  employing  the  first 


moment  PHD  approximation. 

Future  work  will  focus  on  extending  this  solution  to  a  more 
realistic  SSA  scenarios  with  an  arbitrary  number  of  objects 
and  clutter  sources,  both  modeled  using  spatial  Poisson  distri¬ 
butions.  We  also  plan  on  applying  higher  fidelity  six  degree 
of  freedom  dynamics,  modeling  both  three-dimensional  trans¬ 
lation  and  full  rotational  dynamics  as  well  as  (conservative 
and  nonconservative)  perturbations  to  the  two-body  Keplerian 
motion.  Finally,  we  also  seek  to  include  characterization  to  the 
approach  such  that  we  have  a  complete  integrated  detection, 
characterization  and  tracking  AEGIS-FISST  solution  to  the 
general  SSA  problem. 
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